%% Practical: Sensorspace Continuous OAT
%
% In this practical we will run through an analysis step by step, several
% of which require manual intervention. This will go through the following
% steps:
%
%   1) Prepare raw data for OAT analysis
%   2) Bandpass filter data and define a single experimental epoch
%   3) Compute a first level GLM analysis with OAT
%   4) Visualise results with FieldTrip
%   5) Compute a topoplot of a condition contrast across the experiment
%
% This practical we will work with a single subject's data from a blocked
% finger tapping experiment (CTF data courtesy of Matt Brookes), and
% perform a time-frequency analysis in the beta band using a time-wise GLM
% in sensor space.   
% 
% Data can be downloaded from: 
% www.fmrib.ox.ac.uk/~woolrich/ctf_fingertap_subject1_data.tar.gz

%% SETUP THE MATLAB PATHS
%
% Sets the Matlab paths to include OSL. Change these paths so that they 
% correspond to the setup on your computer. You will also need to ensure 
% that fieldtrip and spm are not in your matlab path (as they are included
% within OSL).%% SETUP THE MATLAB PATHS

global OSLDIR;
    
% set this to where you have downloaded OSL:
osldir = '/Users/andrew/Software/Matlab/osl/osl_full/osl1.5.0_beta';  

addpath(osldir);
osl_startup(osldir);

%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS
%
% Change the workingdir variable to correspond to the correct directory.

workingdir='/Users/andrew/Projects/OSL_test/ctf_fingertap_subject1_data'; % directory where the data is

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

cd(workingdir);

%% SET UP THE LIST OF SUBJECTS FOR THE ANALYSIS
%
% Specify a list of the fif files, structual files (not applicable for this
% practical) and SPM files (which will be created). It is important to make
% sure that the order of these lists is consistent across sessions. 
% Note that here we only have 1 subject, but more generally there would be
% more than one, e.g.:

% fif_files{1}=[testdir '/fifs/sub1_face_sss.fif']; 
% fif_files{2}=[testdir '/fifs/sub2_face_sss.fif']; 
% etc...

% spm_files{1} = [workingdir '/sub1_face_sss.mat'];
% spm_files{2} = [workingdir '/sub2_face_sss.mat'];
% etc...
% set up a list of SPM MEEG object file names (we only have one here)
% set up a list of SPM MEEG object file names (we only have one here)

% set up a list of SPM MEEG object file names (we only have one here)
spm_files={[workingdir '/dsubject1.mat']};
 
cleanup_files=0; % flag to indicate that you want to clean up files that are no longer needed

%% LOAD RAW DATA AND INSPECT SPM OBJECT
% 
% The data is in an SPM format so we can use the spm_eeg_load function to load it in.
% By typing "D" we can inspect its contents, the critical bits to check now are that
% we have a continuous (as opposed to epoched) dataset with 216000 samples at 150Hz.

subnum = 1;             
D = spm_eeg_load(spm_files{subnum});

D

%% INSPECT DATA IN OSLVIEW
%
% Note that this is CTF data and so there is only one channel type (axial
% gradiometers).
% You will see that there are no obvious artefacts to mark.

oslview(D);


%% SETUP SENSOR SPACE OAT SOURCE RECON
% 
% This stage sets up the source reconstruction stage of an OAT analysis.
% The source_recon stage is always run even for a sensorspace analysis,
% though in these cases it simply prepares the data for subsequent
% analysis.
%
% In this example we define our input files (only D_continuous) 
% and conditions before setting a time frequency window from  3 to 30Hz and
% 300 to 960 seconds. This time window is the length of the entire task. 
% The source recon method is set to 'none' as we are performing a sensorspace analysis

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Undefined'};
oat.source_recon.freq_range=[13 30]; % frequency range in Hz
oat.source_recon.time_range=[300,32*30];
oat.source_recon.method='none';
oat.source_recon.modalities={'MEGGRAD'};
oat.source_recon.work_in_pca_subspace=1;
oat.source_recon.pca_dim=270;

oat.source_recon.dirname=[workingdir '/subj1_results.oat'];

oat = osl_check_oat(oat);

%% RUN SOURCE RECON
%
% Note that this is a sensor space analysis and so this just preps the data
% for the first level stage  

oat.to_do=[1 0 0 0];

oat = osl_run_oat(oat);

%% ESTABLISH REGRESSOR FOR CONTINUOUS TIME GLM.
%
% This section generates the predictors for our GLM analysis based on the order
% of the finger tapping blocks. There were four conditions:
% Left hand tap
% Right hand tap
% Rest
% Both hands tap
% Rest at start
%
% The order of these blocks in recorded in block_order and their length in seconds
% in block_length.
%
% This should be setup to correspond to the time window
% over which the source recon is carried out.

D=spm_eeg_load(spm_files{1});

x=zeros(length(D.time),5);

block_length=30; %s
block_order=[5 5 5 5 5 5 5 5 5 5 4 3 2 1 2 3 1 4 3 4 1 3 2 1 4 4 2 1 3 3 4 1 4 3 1 2 1 2 3 4 3 4 1 2 3 4 1 2];

% [Left, Right, Rest, Both, Rest_at_start]
% [  1     2      4     8     16 ]
% figure;plot(D.time,squeeze(D(1,:,:)))
% emacs ~/vols_data/From_Nottingham_with_Love/JRH_MotorCon_20100429_01_FORMARK.ds/MarkerFile.mrk 

tres=1/(D.fsample);
tim=1;
for tt=1:length(block_order),    
    x(tim:tim+block_length/tres-1,block_order(tt))=1;
    tim=tim+block_length/tres;
end;

figure;imagesc([1:5],D.time,x);
plot4paper('Regressor','time (secs)');title('Design Matrix');colorbar;

%% SET UP FIRST LEVEL GLM ANALYSIS CONTRASTS
%
% This cell defines the contrast vectors for the first level analysis.
%
% Each contrast is a vector containing a weight per condition defining how
% the condition parameter estimates are to be compared. Each vector will
% produce a different t-map across the sensors. Contrasts 1 and 2 describe
% positive correlations between each sensor's activity and the presence
% of a motorbike or face stimulus respectively. Contrast 3 tests whether
% each sensors activity is larger for faces than motorbikes.
%
% Several other parameters are set afterwards defining the time-frequency
% methods to be used in the first-level analysis, downsampling factor and
% the name of this first-level analysis.
% Note, to look for changes in power specifically in the beta band 
% (since the motor system is known to produce beta band changes) 
% we have set "oat.first_level.tf_freq_range=[13 30];"

oat=osl_load_oat(oat.source_recon.dirname);

oat.first_level.time_moving_av_win_size=1;
oat.first_level.design_matrix=x';
oat.first_level.contrast{1}=[-1 0 1 0 0]'; % rest-left
oat.first_level.contrast{2}=[0 -1 1 0 0]'; % rest-right
oat.first_level.contrast{3}=[0  0 1 -1 0]'; % rest-both
oat.first_level.contrast_name{1}='rest-left';
oat.first_level.contrast_name{2}='rest-right';
oat.first_level.contrast_name{3}='rest-both';
oat.first_level.bc=[0 0 0];
oat.first_level.results.first_level_cons_to_do=1:3;

oat.first_level.doGLM=1;

% Time frequency parameters
oat.first_level.tf_freq_range=[13 30];
oat.first_level.tf_hilbert_freq_res=diff(oat.first_level.tf_freq_range);
oat.first_level.tf_method='hilbert';

oat.first_level.post_tf_downsample_factor=2; 

oat.first_level.name=['subj1_first_level_ft'];

%% RUN FIRST-LEVEL CONTINUOUS TIME OAT ANALYSIS
%
% Run OAT analysis.

oat.to_do=[0 1 0 0];

oat = osl_check_oat(oat);

oat = osl_run_oat(oat);

%% VISUALISE USING FIELDTRIP
%
% Next we will use an osl wrapper around a Fieldtrip function to visualise 
% the results from contrast 3 (Rest>Left) - i.e.,"where in sensor space is 
% the increased beta power in rest compared with left finger tapping?". This 
% requires us to define several parameters in the S2 structure. Critically, 
% we define the oat struct and the index of the first level contrast of interest within it.

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.modality='MEGGRAD';
S2.first_level_contrast=1;
S2.cfg.colorbar='yes';
S2.cfg.interactive='no';
S2.view_cope=0; % means we view the t-stats ( COPE / VARCOPE )

% Talculate t-stat using contrast of absolute value of parameter estimates
osl_stats_multiplotER(S2);

% Try taking a look at the other contrasts and look at the lateralisation


%% LOOK AT BETA POWER AND REGRESSORS AT A SINGLE SENSOR
%
% This section plots the Hilbert envelope of the beta filtered time series
% over time alongside the regressors. The raw time series is extracted by
% rerunning the first-level analysis with the oat.first_level.doGLM flag
% set to zero (returning the time series the GLM would have been fitted
% on).
%
% The sensor with the maximum value in the COPE (Comparison Of Parameter
% Estimates) from the previous GLM analysis.

% Load previous results and select channel of interest
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
[a, chan_ind]=max(squeeze(stats.cope(:,1,S2.first_level_contrast)));
disp(D.chanlabels(chan_ind));

% re-run oat without doing GLM - this will output the time series used to
% fit the GLM to, i.e. the beta power (Hilbert envelope) time courses
oat.to_do=[0 1 0 0];
oat.first_level.doGLM=0;
oat = osl_run_oat(oat);

% Create plot
stats2=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
figure;plot(stats2.glm_input_times,normalise(squeeze(stats2.glm_input_data(chan_ind,:))));

% Compare to design matrix:
time_ind=intersect(find(D.time>=oat.source_recon.time_range(1)),find(D.time<=oat.source_recon.time_range(2)));
ho;plot(D.time(time_ind),x(time_ind,1),'r','LineWidth',2);
plot(D.time(time_ind),x(time_ind,2),'g','LineWidth',2);
plot(D.time(time_ind),x(time_ind,4),'k','LineWidth',2);
legend('data','left','right','both');
plot4paper('time(secs)','beta power');
